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Abstract. - We present a theory for a superfluid Fermi gas near the BCS-BEC crossover, 
including pairing fluctuation contributions to the free energy similar to that considered by 
Nozieres and Schmitt-Rink for the normal phase. In the strong coupling limit, our theory 
is able to recover the Bogoliubov theory of a weakly interacting Bose gas with a molecular 
scattering length very close to the known exact result. We compare our results with recent 
Quantum Monte Carlo simulations both for the ground state and at finite temperature. Excel- 
lent agreement is found for all interaction strengths where simulation results are available. 



The recent experimental realization of strongly interacting Fermi gases of 6 Li and 40 K 
atoms near a Feshbach resonance has opened up the exciting possibility of investigating the 
crossover from a Bardeen-Cooper-Schrieffer (BCS) superfluid to a Bose-Einstein condensate 
(BEC) [1-5]. In these systems, the inter-atomic interaction strength can be varied by tuning 
the energy of a near-resonant molecular state with a magnetic field. 

Below resonance where the s-wave scattering length a is positive, stable diatomic molecules 
are observed to form a BEC at low temperatures. Above resonance, with a < 0, the molecules 
dissociate and form a BCS superfluid of fermionic pairs. In the crossover region where the 
scattering length a is large one can access a new, strongly correlated regime known as the 
unitary limit [6]. Recent experiments in the crossover regime have found evidence for this 
transition by measuring low-lying collective modes [2,7,8] and heat capacity [4,9]. 

These rapid experimental developments constitute an ideal testing ground for theoretical 
studies of the BCS-BEC crossover. However, theoretical results available in the literature are 
limited in the strongly correlated unitary regime. The first systematic study of the crossover 
at zero temperature was provided by Eagles and Leggett based on BCS mean-field equations 
[10, 11]. Later, the effects of pair fluctuations were considered by Nozieres and Schmitt-Rink 
(NSR) at temperatures above the superfluid transition [12,13]. This was recently extended to 
the superfluid phase by Strinati et al. using finite temperature Green functions [14,15]. 

Extensions of these approaches to take into account the bare Feshbach molecule have 
also been presented [16, 17], with the conclusion that additional two-channel effects can be 
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neglected for broad resonances and detunings near the crossover regime [18,19]. All these 
studies give a qualitative description of the crossover: but none of them are quantitatively 
correct in the unitary limit, and in the BEC region. A major drawback of these theories is 
that the predicted value of the molecular scattering length in the deep BEC regime, a m = 2a, 
does not agree with the exact result from the solution of the four-body problem [20-22], 
i.e., a m ~ 0.60a. This much lower value is, however, consistent with Quantum Monte Carlo 
(QMC) simulations [23,24]. 

The purpose of the present Letter is to develop a quantitatively reliable theory for supcrfluid 
Fermi gases in the broad resonance or single channel limit, at all interaction strengths and 
low enough temperatures. To this end, we extend the NSR analysis to the superfluid phase 
on top of the BCS mean field approximation. An essential ingredient of our theory is that 
the number equation, i.e., the relation n = —dCl/d^i, is satisfied for the full thermodynamic 
potential il — not just for the mean-field contribution to 0. In the deep BEC limit, where 
molecule-molecule correlations are important, this requirement renormalizes the mean-field 
molecular scattering length of a m = 2a to a value of a m ~ 0.57a [25], which is very close to 
the exact four-body prediction [20]. As a consequence, our results for the equation of state at 
zero temperature along the full range of the BCS-BEC crossover are in excellent agreement 
with the QMC data [23]. The temperature dependent results also agree with recent path 
integral Monte Carlo calculations [26]. 

The system we consider is a uniform gas of N Fermi atoms in two hyperfine states denoted 
as pseudo-spins a =\, \, with = N± = N/2. To characterize the superfluid ground state, 
we introduce explicitly an order parameter A that will be determined at the mean field level, 
and use the Nambu spinor representation, in which the system is described by the Hamiltonian 
H = Ho + Va + Vm, with the terms, 

n ° = X^k [&<7*-ACT x ]tfk + 2&, (!) 

k k 

Va = ^^A<7 x Vk, (2) 

k 

u 

kfe'q 

where ip£ = (c^f, c_ki) is the Nambu creation field operator for Fermi atoms with the kinetic 
energy £k = £k — A* = ?i 2 k 2 /2m — fi, fj, is the chemical potential, and a x and a z are the 2x2 
Pauli matrices. The contact interaction Uo for atom-atom interactions is renormalized by 
introducing the s-wave scattering length a. This is defined by the low energy limit of the 
two-body scattering problem via (47r7i 2 a/m) _1 = U ~ 1 + X)k(2 £ k) _1 - 

After taking into account anti-commutators, the "unperturbed" term Ho can be identified 
as the standard BCS mean field Hamiltonian. In variational calculations, the approximate 
ground state would be an eigenstate of Ho, with A varied to minimize the total energy [16]. 
However, this approach is inappropriate on the BEC side of the transition, where molecule- 
molecule scattering gives rise to strong density correlations involving four particles rather 
than just two. The purpose of the present Letter is to describe an approximate perturbation 
expansion which takes this into account. 

A diagrammatic expansion is performed in terms of the propagator of the "free" Hamilto- 
nian Ho, which is given by: 



Vint = Y Yl (Vk +9 ^k) (Vi-c^k') • (3) 
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Fig. 1 - Diagrammatic representations of the mean field contributions [(a) and (b)] and the pairing 
fluctuation contributions (c) to the thermodynamic potential. The full line represents Go (k, iu m ), 
the line with arrows describes the anomalous Green function G12 = G21 = J- (k,iui m ), and the dotted 
line is Uo- 



Here Q{k,iu m ) = -[iu m + { k ]/[^ + E%\ and T(\t,iu} m ) = A/[cj^ + E%\ are, respectively, 
the BCS normal and anomalous Green functions, while = + A 2 ] 1 / 2 denotes the single 
particle excitation energy and u> m — (2m+l)7r/(7i/3) is the fcrmionic Matsubara frequency with 
inverse temperature (3 = 1/ksT. The associated "unperturbed" thermodynamic potential 
then takes the form, fi„ = l//3£ km TrlnG (k, iw m )+£ k & - E k [(£k-£k)+2//31n /(-£ k )], 
where f(x) = \j(eP x + 1) is the Fermi function. 

The contributions of interactions to the thermodynamic potential consist of a static mean 
field part and a fluctuation part originating from the particle-particle Cooper channel. As 
shown diagrammatically in Fig. la and Fig. lb, the mean-field corrections from Va and Vj n t 
read, respectively, Sil A = l//3£k,m ATr[a :E G (k, itu m )} = -2A 2 /U and 8tt mf = U [l//3 
Ek m 3~0*-, iu m )] 2 = A 2 /{/q, where we have used the mean field gap equation: 
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Here the coupling Uo is eliminated in favor of the s-wave scattering length a. These corrections, 
together with f^o, give rise to an overall mean field thermodynamic potential of £l m f = + 

^/=Et^-^k + ^ + |ln/(- J B k )]--P^. (6) 
k ^e k p inn a 

Fluctuation corrections beyond mean field are illustrated in Fig. lc. A sum of the resulting 
geometrical series thus leads to [12, 13, 17], 



1 f +oa 1 

°" f = --E/ ^^T^r^q^)' ( 7 ) 



7T 

q 



where, following NSR, we have written fl p f in terms of a phase shift defined by <5(q, ui) = 
- Im ln[-xn(q, a + irjj\ - 1/2 Imln{l - xi 2 (q> ^ + "?)/[Xn(q, <*> + "))Xii(q, -w + «??)]}. Here 
the analytic continuation is performed and 77 is a positive infinitesimal, while xn = l/^o + 
V/'Ek.m^l" fc,zz/„ - iu; m )5(k,iw m ) and X12 = l//3£)k,m -^(q - - iw m )J"(k, 

are the diagonal and off-diagonal parts of the Cooper pair propagator, with v n — 2mr/P being 
the bosonic Matsubara frequency. 

Putting together the mean field and the pairing fluctuation corrections to the thermo- 
dynamic potential, we obtain the total contributions, Q = Q m f + £l p f. We emphasize that 
the gap (J5J is chosen at the mean field level, i.e., d£l m f/dA = 0. Non-trivial effects be- 
yond the BCS mean field approximation enter into the theory through the modified num- 
ber equation N = — (dSl/d[i) T . Explicitly, we obtain N = N m f + N p f ttJj + N p f t A, where 
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N m f = -(dn m //d/z)rA, Npf.f, = -(dfl p f/d/j,) T A, and A^/.a = -(d£l P f/dA) Tli (dA/dfi). 
The coupled gap and particle number equations, together with the thermodynamic potentials 
© and Q, form the basis of our NSR theory in the broken-symmetry state. Once \x and 
A are obtained as a function of the interaction strength and temperature, the entropy and 
the energy of the gas can then be calculated straightforwardly, using S = — (dCl/dT) and 
E = Q + TS + fj,N. 

There is a key difference between our method and the diagrammatic theory proposed by 
Strinati et al. [14,15]. The term N p f^, which ensures number conservation, is not included 
in the diagrams considered by Strinati et al. This term becomes increasingly important in 
the BEC regime. Physically, it generates the dominant part of the four-fermion correlations, 
which are increasingly important for tightly bound Cooper pairs. This constitutes a major 
advantage of our NSR theory. 

To give more insight, in the strong coupling limit we re-interpret our formalism in the 
framework of the functional integral method [13], in which a Cooper-pair Bose field A(x, r) is 
introduced through the Hubbard-Stratonovich transformation. Integrating out the fermionic 
degrees of freedom in the usual fashion and setting A(x, r) = A + 8A(x, r), the resulting ef- 
fective bosonic action is then expanded up to quadratic order in fluctuations 8A(x, r): S e ft rs 
SW+S&h After performing a Fourier transformation we find that = [3A 2 /Uq— Trln[— G^" 1 ]-!- 
/?Ek6c, and 

Si2) = EM^Xiifa ii/„)*A, + X i2(q, ip n )(5A+SA± q + 5A_ q 5A q )/2] , (8) 

q.n 

where 8A q = <5A(q, iv n ) is the Fourier transformation of 5A(x,t). The saddle point so- 
lution of gives the standard mean-field theory, while the next order Gaussian expan- 
sion leads to exactly the same contribution as in Eq. J2J). It is easy to see that in 
the long-wavelength and low-frequency limits, — Xn(q, iv n ) ~ —iv n + ?i 2 q 2 /4m + A 2 / (— 4/z) 
and Xi2(q, ivn) ~ A 2 / (— 4/x). Hence, acquires the familiar form of Bogoliubov excita- 
tions [27]. Our formalism therefore incorporates interactions between condensed and non- 
condensed Cooper-pairs at the level of Bogoliubov theory, which must be the dominant part 
of four-fermion correlations at low temperature. 

The crucial observation of the present Letter is that in the deep BEC limit (A/(— fi) — > 0) 
N p f t A ~ A 2 /(— /i) 1 / 2 is of the same order of N m f, or more precisely, N p f t A = CN m f with 
C ~ 2.5, as shown analytically in a forthcoming publication [28]. In contrast, N p f tfl ~ 
A 3 /(— [ifl" 1 ~ A/(—fj,)N m f becomes negligible. We thereby find that the molecular conden- 
sate density Ng ~ A^/2 ~ (1 + C)N m f/2. Mean-field number and gap equations provide the 
expressions, N m f ~ A 2 m 2 a/(47r?i 2 ), and fig = 2/i + h 2 /(ma 2 ) ~ A 2 ma 2 / (2h 2 ), where /ie is 
the molecular chemical potential. Assembling together these three expressions and eliminating 
A, one finds that [i B = 4ttH 2 [2o/(1 + C)] / (2m) x N%, implying a m = 2a/ (I + C) ~ 0.57a - 
a value close to the exact result (a m ~ 0.60a). The residual 5% difference occurs because we 
exclude interactions between non-condensed Cooper-pairs. 

Our diagrammatic procedure is also different from the works in Refs. [20, 22], where 
molecule-molecule scattering processes are considered in isolation. These calculations per- 
turb around the free-fermion state and therefore require additional terms beyond the ladder 
structures in Fig. 1 to obtain correct results. By comparison, we perturb around a BCS 
state of correlated fermions below T c . A dominant part of the molecule- molecule scattering is 
therefore already included in the ladder diagrams due to the inclusion of N p f^, as we have 
shown using the functional integral method. To contain the full four-body correlations, a more 
complicated wave function has to be implemented [21]. Equivalently, within the functional 



H. Hu, X.-J Liu and P. D. Drummond: EQUATION OF STATE IN THE BCS-BEC CROSSOVER!) 




Fig. 2 - (color online) Energy per particle E/N, in the BCS-BEC crossover (with the binding energy 
65 = —ft 2 /ma 2 subtracted). Lines plotted are: our result (solid line); fixed node QMC data from 
Ref. [23] (circles with error bars), Lee and Yang perturbation theory (dashed line); BCS mean field 
predictions (dot-dashed line). Inset: enlarged view of the BEC regime. The solid line corresponds 
to the equation of state of a repulsive gas of molecules given by LHY [30], {E/N — Eb/2)/eFG = 
5/(187r)fc F am[l + 128/(15v / 67r3)(fc F (i m ) 3/2 + ...], with our analytic result of a m = 0.57a. 



integral method, we may include four-fermion correlations by expanding the action to fourth 
order of the gap [25]. However, this calculation is difficult to extend to the crossover regime. 

Figure 2 presents our results for the energy per particle as a function of the interaction 
strength 1/fcpci at a low temperature T = 0.02Tp. The energy scale is given by the noninteract- 
ing energy, 6fg = (3/5)ep = (3/lO)h 2 k F /m, where Uf = (S^n) 1 / 3 is the Fermi wave vector. 
As a benchmark, the approximate fixed node QMC data at zero temperature is shown [23] , to- 




Fig. 3 - (color online) Chemical potential in the BCS-BEC crossover, as predicted by various ap- 
proaches. The QMC curve in dotted line is calculated from a best fit to the QMC energies, as 
outlined in Ref. [23] . The diagrammatic predictions (empty squares) by Strinati et al. are taken from 
Ref. [15] without the inclusion of the self-energy shift. 
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Fig. 4 - (color online) Temperature dependence of the energy per particle (main plot) and the chemical 
potential (inset) in the unitary limit. Lines plotted are: our results (solid lines); path integral Monte 
Carlo data from Ref. [26] (circles); BCS mean field (dot-dashed lines); the original NSR prediction in 
Ref. [12,131. 



gether with the perturbation results E/(Ne FG ) = l + 10/(97r)A: F a+4/(2l7r 2 )(ll-2 ln2)(fc F a) 2 
obtained by Lee and Yang [29]. In the BCS region, 1/kpa < —0.5, we find that our results 
agree with the QMC data, apart from residual differences due to the energy dependence of the 
scattering amplitude. In the BEC region, 1/kpa > 0.5, our predictions coincide perfectly with 
QMC calculations, and thereby agree with the equation of state of a repulsive gas of molecules 

with a m ~ 0.60a, as shown in the inset. In the range -0.5 ~ \/{k F a) ~ 0.5 spanning the 
most interesting crossover region, our results differ only slightly from that of QMC simula- 
tions. In particular, in the unitary limit we predict E/N — £cfg with £ = 0.401, compatible 
with the QMC findings £ = 0.42(1) [23]. The overall agreement between the two alternative 
calculations is therefore excellent, especially in the challenging strong coupling regime. 

We have also calculated the chemical potential at the BCS-BEC crossover, predicted by 
different theories (Fig. 3). We find again an excellent quantitative agreement between our 
NSR results and QMC calculations. To emphasize the similarity between our formalism and 
the diagrammatic theory given by Strinati et al. [14], we compared the results obtained without 
the inclusion of the crucial term N p f a in the number equation, with Strinati's diagrammatic 
findings [15]. These asymptotically approach the mean field predictions in the deep BEC 
limit, which would (incorrectly) imply that a m = 2a. This observation thereby unambiguously 
verifies that N p f a is responsible for obtaining the correct molecular scattering length. 

Finally, in Fig. 4, we compare our predictions at finite temperature with recent path inte- 
gral Monte Carlo calculations of spin 1/2 fermions [26] in the strongly coupled unitary limit. 
At low temperature up to T c « 0.22TV, these results are in good qualitative agreement with 
each other. The residual discrepancy is possibly due to finite size effects in the simulations. 

To conclude, we have presented an NSR-type formalism for a Fermi gas at the BCS- 
BEC crossover, in the broken-symmetry phase. A notable achievement of our formalism is 
that a Bogoliubov theory of composite Cooper pairs is reproduced in the BEC limit, with a 
molecular scattering length very close to the exact value. We have compared our predictions 
of the equation of state of the gas with available Monte Carlo calculations, and find excellent 
agreement. Our results also make quantitative contact with a previous diagrammatic theory in 
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the weak and intermediate coupling regimes. We believe, therefore, that the present formalism 
provide a quantitatively reliable description of superfluid Fermi gases at low temperature over 
the entire range of the crossover. 
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